clc
clear

z_max = 40000;
nc    = 40; 
% number of cells on vertical distribution
m     = 2.5;    % stretch coefficient on x-dir

nl    = nc+1; % number of levels(cell boundaries) on vertical distribution
eta_max = nc;
deta  = 1;
eta   = 0:deta:eta_max;
eta   = eta';

x = m * pi * 2 * ( eta ./ eta_max - 0.5 );
k = z_max * 2; %  stretch coefficient on y-dir
f = k * ( atan(x) + pi/2 ) / pi;
df = k ./ ( 1 + x.^2 )/ pi;

F_n = zeros(nl,1);
for i = 2:nl
    F_n(i) = F_n(i-1) + f(i) * deta / eta_max;
end

x_min = -pi * m;
x_max = x;
F_a = k / ( 2 * pi^2 * m ) * ( ( x_max .* atan(x_max) - log( sqrt( 1 + x_max.^2 ) ) ) - ( x_min .* atan(x_min) - log( sqrt( 1 + x_min.^2 ) ) ) ) + k / 2 * eta ./ eta_max ;

figure
hold on
plot(eta,f,'r')
plot(eta,f,'Color','r','Marker','o')
plot(eta,F_a,'b')
plot(eta,F_a,'Color','b','Marker','o')
plot(eta,df,'k')
plot(eta,df,'Color','k','Marker','o')

level = zeros(2,nl);
for i = 1:nl
    level(:,i) = F_a(i);
end

figure
plot(level,'b')